package fem2.material;

import inf.jlinalg.ArrayVector;
import inf.jlinalg.BLAM;
import inf.jlinalg.IMatrix;
import inf.jlinalg.IVector;

import java.util.ArrayList;

import math2.LinAlgUtilities;
import math2.MathUtilities;
import math2.MutableDouble;
import fem2.Element;
import fem2.MeshPart;
import fem2.NotImplementedError;
import fem2.Variable;
import fem2.element.StructuralElement;
import fem2.enu.State;
import fem2.enu.VariableName;

/**
 * Simplified Lemaitre elastoplastic damage material model
 * 
 * <pre>
 * Reference: Souza Netp, Computational Plasticity, section 12.4
 * </pre>
 * 
 * TODO: need check
 * 
 * @author hbui
 * 
 */
public class LemaitreDamageMM2 extends InelasticMaterialModel {

	protected FluencyCriteria isotropicLaw;
	protected double r, s;// material constants

	public LemaitreDamageMM2(double E, double nu, double thickness, double rho, State ss,
			FluencyCriteria isotropicLaw, double r, double s) {
		super(E, nu, thickness, rho, ss);

		if (getState() == State.PLANE_STRESS) {
			throw new Error(
					"Because of the complexity, this material model does not support plane stress version. Please use the plane stress version instead");
		}

		this.isotropicLaw = isotropicLaw;
		this.r = r;
		this.s = s;
	}

	@Override
	public void Initialize(Element e) {
		super.Initialize(e);

		int dim = e.getMeshPart().getDimension();
		int ntrs = dim * (dim + 1) / 2;
		int ntrs1 = (dim == 2 ? ntrs + 1 : ntrs);

		ArrayList<Variable> vars = getVariables(e);
		for (Variable v : vars) {
			v.registerScalarVariable(VariableName.ALPHA, true);// strain-like
																// variable (R)
			v.registerScalarVariable(VariableName.DAMAGE, true);// damage
																// variable (D)
			v.registerVectorVariable(VariableName.PLASTIC_STRAIN, ntrs1, true);
			/*
			 * save the variable to create history
			 */
			v.save();
		}
	}

	@Override
	public void InitalizeNonlinearIteration(Element e) {
	}

	@Override
	public void FinalizeNonlinearIteration(Element e) {
		StructuralElement se = (StructuralElement) e;
		MeshPart mp = e.getMeshPart();
		int dim = mp.getDimension();
		int ntrs = dim * (dim + 1) / 2;
		int ntrs1 = (dim == 2 ? ntrs + 1 : ntrs);
		int ng = e.getNumIntegrationPoints();

		IVector eetrial = new ArrayVector(ntrs1);
		IVector eedtrial = new ArrayVector(ntrs1);
		IVector s = new ArrayVector(ntrs1);

		if (dim == 2) {

			/*************************************************************
			 * Von Mises plane strain
			 *************************************************************/

			for (int k = 0; k < ng; k++) {
				Variable var = getVariables(e).get(k);
				double[] xi = e.getIntegrationPoint(k);
				IVector sigma = se.getStress(k);
				IMatrix C = getConstitutiveMatrix(e, k);
				IVector epsilon = var.getVectorVariable(VariableName.STRAIN);
				IVector ep = var.getVectorVariable(VariableName.PLASTIC_STRAIN);
				IVector ep_old = var.getVectorVariable(VariableName.PLASTIC_STRAIN, -1);
				MutableDouble damage = var.getScalarVariable(VariableName.DAMAGE);
				double damage_old = var.getScalarVariable(VariableName.DAMAGE, -1);
				MutableDouble R = var.getScalarVariable(VariableName.ALPHA);
				double R_old = var.getScalarVariable(VariableName.ALPHA, -1);
				double w_old = 1 - damage_old;

				/*
				 * compute trial elastic strain
				 */
				BLAM.copy(se.computeStrain(xi), epsilon);
				eetrial.set(0, epsilon.get(0) - ep_old.get(0));
				eetrial.set(1, epsilon.get(1) - ep_old.get(1));
				eetrial.set(2, epsilon.get(2) - ep_old.get(2));
				eetrial.set(3, -ep_old.get(3));

				/*
				 * compute trial volumetric elastic strain
				 */
				double eevtrial = eetrial.get(0) + eetrial.get(1) + eetrial.get(3);

				/*
				 * compute trial deviatoric elastic strain
				 */
				eedtrial.set(0, eetrial.get(0) - eevtrial / 3);
				eedtrial.set(1, eetrial.get(1) - eevtrial / 3);
				eedtrial.set(2, eetrial.get(2) / 2); // physical value
				eedtrial.set(3, eetrial.get(3) - eevtrial / 3);

				/*
				 * compute trial effective hydrostatic stress
				 */
				double ptilde = kappa * eevtrial;

				/*
				 * compute equivalent stress
				 */
				double J2 = Math.pow(2 * mu, 2)
						* (Math.pow(eedtrial.get(2), 2) + 0.5 * (Math.pow(eedtrial.get(0), 2)
								+ Math.pow(eedtrial.get(1), 2) + Math.pow(eedtrial.get(3), 2)));
				double qttrial = Math.sqrt(3 * J2);

				/*
				 * compute yield stress
				 */
				double oY = isotropicLaw.computeFluencyValue(R_old);

				/*
				 * check plastic admissibility
				 */
				double phi = qttrial - oY;
				if (phi <= 0.0) {

					BLAM.copy(ep_old, ep);

					damage.setValue(damage_old);

					R.setValue(R_old);

					double w = 1 - damage.getValue();

					sigma.set(0, w * (2 * mu * eedtrial.get(0) + ptilde));
					sigma.set(1, w * (2 * mu * eedtrial.get(1) + ptilde));
					sigma.set(2, w * 2 * mu * eedtrial.get(2));

					BLAM.copy(De, C);
					BLAM.scale(w, C);

				} else {

					/*
					 * compute delta lambda
					 */
					double dlambda0 = phi * w_old / (3 * mu);
					double dlambda = computeDeltaLambda(dlambda0, ptilde, qttrial, R_old, w_old);

					/*
					 * update internal variable
					 */
					R.setValue(R_old + dlambda);

					/*
					 * compute new yield stress
					 */
					oY = isotropicLaw.computeFluencyValue(R.getValue());

					/*
					 * compute material integrity (eq 12.46)
					 */
					double w = 3 * mu * dlambda / (qttrial - oY);

					/*
					 * update equivalent stress
					 */
					double q = w * oY;

					/*
					 * update hydrostatic stress.
					 */
					double p = w * ptilde;

					/*
					 * update deviatoric stress.
					 */
					double factor = 2 * mu * q / qttrial;
					s.set(0, factor * eedtrial.get(0));
					s.set(1, factor * eedtrial.get(1));
					s.set(2, factor * eedtrial.get(2));
					s.set(3, factor * eedtrial.get(3));

					/*
					 * update sigma
					 */
					sigma.set(0, s.get(0) + p);
					sigma.set(1, s.get(1) + p);
					sigma.set(2, s.get(2));

					/*
					 * update plastic strain
					 */
					factor = 1 - 3 * mu * dlambda / (w * qttrial);
					eetrial.set(0, factor * eedtrial.get(0) + eevtrial / 3);
					eetrial.set(1, factor * eedtrial.get(1) + eevtrial / 3);
					eetrial.set(2, factor * eedtrial.get(2) * 2);
					eetrial.set(3, factor * eedtrial.get(3) + eevtrial / 3);

					ep.set(0, epsilon.get(0) - eetrial.get(0));
					ep.set(1, epsilon.get(1) - eetrial.get(1));
					ep.set(2, epsilon.get(2) - eetrial.get(2));
					ep.set(3, -eetrial.get(3));

					/*
					 * update damage variable
					 */
					damage.setValue(1 - w);

					/*
					 * compute thermodynamics force (eq 12.47)
					 */
					double Y = -(Math.pow(oY, 2) / (6 * mu) + Math.pow(ptilde, 2) / (2 * kappa));

					/*
					 * update tangent operator
					 */
					double H = isotropicLaw.computeFluencyDerivative(R.getValue());
					double wp = (3 * mu + w * H) / (qttrial - oY);
					double Fp = wp - (H / (3 * mu)) * Math.pow(-Y / r, this.s)
							+ (this.s * H * oY * (qttrial - oY) / (9 * Math.pow(mu, 2) * r))
							* Math.pow(-Y / r, this.s - 1);

					double a1 = (w / (qttrial - oY) - Math.pow(-Y / r, this.s) / (3 * mu)) / Fp;
					double a2 = -(this.s * ptilde * (qttrial - oY) / (3 * mu * r * kappa * Fp))
							* Math.pow(-Y / r, this.s - 1);
					double a3 = a2 * wp;
					double a4 = a1 * wp - w / (qttrial - oY);

					double aFactor = 2 * mu * w * oY / qttrial;
					double bFactor = 2 * mu * (a1 * H * w + a4 * oY - w * oY / qttrial);
					double cFactor = kappa * sqrt_2_3 * (a2 * H * w + a3 * oY);
					double dFactor = 2 * mu * sqrt_3_2 * ptilde * a4;
					double eFactor = kappa * (w + a3 * ptilde);

					double norm_s = Math.sqrt(Math.pow(s.get(0), 2) + Math.pow(s.get(1), 2) + 2
							* Math.pow(s.get(2), 2) + Math.pow(s.get(3), 2));
					BLAM.scale(1 / norm_s, s);

					for (int i = 0; i < ntrs; i++) {
						for (int j = 0; j < ntrs; j++) {
							C.set(i, j,
									aFactor * LinAlgUtilities.Idev2d.get(i, j) + bFactor * s.get(i)
											* s.get(j) + cFactor * s.get(i)
											* LinAlgUtilities.SOID[j] + dFactor
											* LinAlgUtilities.SOID[i] * s.get(j) + eFactor
											* LinAlgUtilities.OneByOne2d.get(i, j));
						}
					}
				}
			}
		}

		else if (dim == 3) {

			/*************************************************************
			 * Von Mises three dimensional
			 *************************************************************/

			// for (int k = 0; k < ng; k++) {
			// double[] xi = e.getGaussPoint(k);
			// }

			/*
			 * TODO
			 */

			throw new NotImplementedError();
		}
	}

	protected double computeDeltaLambda(double dlambda0, double ptilde, double qttrial,
			double R_old, double w_old) {
		double dlambda = dlambda0;
		// Debugger.watch("dlambda0 =", dlambda);
		int cnt = 0;
		while (++cnt < MAX_NT) {
			double oY = isotropicLaw.computeFluencyValue(R_old + dlambda);
			double Y = -(Math.pow(oY, 2) / (6 * mu) + Math.pow(ptilde, 2) / (2 * kappa));
			double tmp = 3 * mu / (qttrial - oY);
			double w = tmp * dlambda;
			double F = w - w_old + Math.pow(-Y / r, s) / tmp;
			// Debugger.watch("F =", F);

			if (Math.abs(F) < MathUtilities.tol6) {
				break;
			}

			double H = isotropicLaw.computeFluencyDerivative(R_old + dlambda);
			double dY = -H * oY / (3 * mu);
			double dF = tmp + tmp * dlambda * H / (qttrial - oY) - H / (3 * mu)
					* Math.pow(-Y / r, s) - (s / r) * (dY / tmp) * Math.pow(-Y / r, s - 1);

			dlambda += -F / dF;

			// Debugger.watch("dlambda =", dlambda);
		}
		// Debugger.watch("----");
		if (cnt == MAX_NT) {
			throw new Error("dlambda computation does not converge in " + MAX_NT + " loops");
		}
		return dlambda;
	}

	// protected double computeDeltaLambda(final double dlambda0, final double
	// ptilde,
	// final double qttrial, final double R_old, final double w_old) {
	//
	// UnivariateFunction function = new UnivariateFunction() {
	//
	// @Override
	// public double value(double x) {
	// double oY = isotropicLaw.computeFluencyValue(R_old + x);
	// double Y = -(Math.pow(oY, 2) / (6 * mu) + Math.pow(ptilde, 2) / (2 *
	// kappa));
	// double tmp = 3 * mu / (qttrial - oY);
	// double w = tmp * x;
	// double F = w - w_old + Math.pow(-Y / r, s) / tmp;
	// return F;
	// }
	// };
	//
	// final double relativeAccuracy = 1.0e-12;
	// final double absoluteAccuracy = 1.0e-8;
	// final int maxOrder = 5;
	// UnivariateSolver solver = new
	// BracketingNthOrderBrentSolver(relativeAccuracy,
	// absoluteAccuracy, maxOrder);
	// double dlambda = solver.solve(100, function, 0, 1e3, dlambda0);
	//
	// Debugger.watch("F = ", function.value(dlambda));
	//
	// return dlambda;
	// }
}
